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Abstract 

The interdiffusion of a solvent into a polymer melt has been studied using large scale molecular 
dynamics and Monte Carlo simulation techniques. The solvent concentration profile and weight 
gain by the polymer have been measured as a function of time. The weight gain is found to scale 
as t^/^, which is expected for Fickian type of diffusion. The concentration profiles are fit very well 
assuming Fick's second law with a constant diffusivity. The diffusivity found from fitting Fick's 
second law is found to be independent of time and equal to the self diffusion constant in the dilute 
solvent limit. We separately calculated the diffusivity as a function of concentration using the 
Darken equation and found that the diffusivity is essentially constant for the concentration range 
relevant for interdiffusion. 
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I. INTRODUCTION 



The interdiffusion of a solvent into a polymer has been a subject of experimental and theo- 
retical research due to both its scientific and practical importance. There has been a number 
of studies of penetrant diffusion in homopolymer a'^i^i^i'^i^i^i'''i^i^i'^°i'^^i^^i^^i^^i^^i^^'^'''i'^^ but less on 
interdiffusion of solvent into polyme r^^i^°i^^i^^i^^ or polymer-polymer interdiffusion.— Pre- 
dicting accurately the nature of the interdiffusion of a solvent into a polymer film has turned 
out to be a challenging problem due to the large number of factors that control the diffusion 
process, including the molecular weight distribution of the polymer and size of the solvent.— 
Whether the polymer is a melt above its glass transition or an amorphous solid below the 
glass transition significantly changes the interdiffusion process. 

For a polymer melt, if a solvent film is placed in contact with the polymer as shown in 
Fig. ^ then the diffusion is one- dimensional and can often be described satisfactorily by 
Pick's second law2^ 

at = fef^Ws*' 

where c is the solvent concentration and D{c) is the diffusivity. This equation assumes 
that the volume of the medium is not changed by the interdiffusion of the solvent. In this 
case the nature of the diffusion is called Pickian or Case I .^^i^^ One fingerprint for Pickian 
diffusion is that the penetration or weight gain by the polymer system increases as t^/^. In 
general D{c) is dependent of concentration c and Eq. Qmust be solved numerically, except 
in special cases.-^ If the diffusivity D{c) = Do is independent of solvent concentration c then 
the solution of Eq. ^ for the concentration of solvent in the medium as a function of time 
and position is simply 

c{z,t) =co{l-eii(^z/2^/lDot)^y (2) 

Here cq is the equilibrium solvent concentration in the polymer usually expressed in units of 
mass per unit volume and erf is the error function. Even though the solvent in general swells 
the polymer and -D(c) may not be independent of c, the simple functional form Eq. (2) is 
often used to fit experimental data.— 

Por glassy polymers the diffusion process does not always follow the standard Pickian 
model and is, in general, referred to be anomalous or non-Pickian diffusion phenomena, which 
is caused by viscoelastic effects in the polymer-solvent system. One type of anomalous behav- 
ior, which has been observed experimentally, is called Case Ij22i22i^2i21i^ in which the polymer 
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FIG. 1: Schematic representation of the MD simulation box for the interdiffusion study. To ahow 
diffusion in one direction the system is periodic in the x and y direction but not in z. The vapor 
region is added to keep the pressure constant. 

relaxation process is very slow compared to the diffusion and exhibits a sharp concentration 
front that propagates at constant speed. However, recent careful experiment3i^>22iSii2^ have 
shown that a Fickian-like precursor foot proceeds the sharp front. To our knowledge no 
simulations have been done to date on Case II diffusion due to the extensive computational 
effort required. In this paper we report on our interdiffusion studies of a solvent into a 
homopolymer above the glass transition temperature, which is expected to exhibit Fickian 
diffusion behavior. We leave the anomalous interdiffusion into a glassy polymer for a later 
study. 

The transport of penetrant molecules in rubbery polymers has been studied using molec- 
ular dynamics (MD) simulation techniques.— i^i'^i^i^i^i^i^^i^'^i'^^i'^'^i'^^i'^^ Most of these studies have 
focused on the penetrant diffusion of solvent molecules in a polymer to determine the self 
diffusion constant. The self diffusion constant -Ds(c) of the solvent is easily calculated from 
the slope of the solvent mean square displacement for long times according to the Einstein 
relation 

. >„„ MM). (3, 

In eqn. El the (...) denote an ensemble average and is obtained by averaging over all solvents 
and many initial time origins. In general -D(c) and -Ds(c) are equal only in the dilute solvent 
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limit, c ~ 0. 

There have been a number of proposed relations between D{c) and Ds{c), the most 
common is the Darken equation^ 

where Dc{c) is called the corrected diffusivity and is related to molecular mobility and / is 
the fugacity of the solvent in the polymer. The thermodynamic factor din f /dine goes to 
unity in the dilute limit. Eq. 0] assumes that diffusion is driven by gradients in chemical 
potential. The corrected diffusivity can be expressed microscopically a o^'^i^^i^^ 

1 

^'^^'^ = 3iV^ Jo ^^^^^ ' 
where J{t) is the interdiffusion current, 

i=l j=l 

Here Xg and Xp are the mole fractions of the solvent and the polymer, respectively and Ng 
and Np are the number of solvent and polymer monomers. The Einstein form of Eq. Eli^^ 

Dc{c) = NxsXp lim ^{{[rcm,s{t) - rcm,s{0)] - [rcm,p(t) - rcr„,p(0)]}^) (7) 

where rcm,s{t) and rcm.,p{t) are the center of mass of all solvent monomers and all polymer 
molecules at time t, respectively. The Einstein form is preferable to the Green-Kubo form 
(velocity auto correlation) since it avoids the need to integrate the time correlation functions. 

To the best of our knowledge, no MD simulation has been carried out to study the 
relationship between the interdiffusion of a solvent into a polymer and the diffusivity and 
self diffusion of solvent in an equilibrated polymer solution. Unlike solvent diffusing in a 
zeolite^i^i^ or other microporous media,-^ the polymer swells as the solvent interpenetrates, 
making it difficult to measure the concentration dependence of the diffusivity D{c) directly. 

Here we are mainly interested in investigating the penetration rates and concentration 
profiles as a function of time through direct analysis of molecular trajectories for Fickian dif- 
fusion in amorphous polymer-solvent systems. We are particularly interested in the relation 
between the diffusivity -D(c), the self diffusion constant Ds{c) and the corrected diffusion 
constant Dc{c) for the solvent. Under these circumstances, MD is a useful tool to deter- 
mine the desired quantities and we have used it for our present study. By combining MD 
simulation with Monte Carlo^ methods we can determine all of these quantities. 
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The paper is organized as follows. In the following section a brief review of the model 
and the simulation method used is given. The two types of thermostats, Langevin and 
Dissipative Particle Dynamics (DPD), used in the simulation are also briefly described. In 
Sec. Ill the results for the diffusion constants as a function of solvent concentration using 
the two thermostats are presented. In Sec. IV the interdiffusion results are presented and 
discussed. Finally, the main results are summarized in Sec. V. 



II. MODEL AND SIMULATION DETAILS 



We performed MD simulations of polymer-solvent system using the coarse grained bead- 
spring model which has been applied successfully to study the effect of entanglement in 
polymer melts.— In this model the polymers are represented by freely jointed bead spring 
chains of length monomers of mass m. The solvent is modeled as either single monomers 
of mass m or dimers of mass 2m. The potential energy associated with interaction between 
nonbonded monomers of type a and (3 is given by the standard Lennard- Jones 6-12 potential 

(8, 

[ 0, r > 

where r is the distance between monomers. Here we take a = 0-^13 and = 2.5a. eaf3 defines 
the units of energy. In this study we set e = epp = ess where p stands for a polymer monomer 
and s for a solvent monomer and vary the relative interaction e^p. 

In addition to the Lennard- Jones interaction between bonded monomers we add an an- 
harmonic interaction term known as FENE potential, 

U^r) = I ~ ^^^^'^'^ ' - ^° (9) 

I 00, r < Rq 

where, as in previous studies^^*^ Rq = 1.5a and k = 30e. 

Our system is also weakly coupled to a heat bath in order to keep it at the desired 
temperature and as a result each monomer moves according to the following stochastic 
equation of motion 



m- 



V$^t/(r.,) + Ff + Ff, (10) 



where m is the monomer mass, U{rij) is the sum of the Lennard- Jones and an harmonic 
spring potential and and Ff are the dissipative force and random force, respectively. 



These later two terms, and together define the type of thermostat used in the 
simulation. In many simulations of polymer melts, a Langevin thermostat is used in which 
the particles are coupled weakly to a heat bath. For polymer melts this is a good way to 
thermostat the system since long range hydrodynamic interactions are screened. However 
these same interactions are important in the diffusion of small molecules as in the pure 
solvent. For this reason we have used two types of thermostats: the Langevin thermostat, 
which by its nature screens the hydrodynamics interactions and the Dissipative Particle 
Dynamics (DPD) thermostat, which does not.-^^^ Part of the aim of this study is to compare 
the results for the two thermostats. 

A. Langevin Thermostat 

In the Langevin thermostat the dissipative part of the force takes the form of friction 
that is proportional to the monomer velocity 

F? ^ -,n,% (U) 

where 7 is the damping constant which is the same for all monomers. The random force is 
related to the frictional force by the fluctuation/dissipation theorem 

(Ff (t) ■ Ff (t')) = ^iksTmS^M - t'y, (Ff (t)) = 0, (12) 

where ks is Boltzmann's constant and T is the temperature. The damping constant 7 con- 
trols both the variance of the random force and the magnitude of the frictional force and 
ensures that the system is kept stable at the desired temperature through out the simula- 
tion. However, since the frictional force does not conserve momentum, the hydrodynamic 
interactions are screened. 

B. Dissipative Particle Dynamics (DPD) Thermostat 

DPD has been first introduced^^ and applied to different system a^^i"^'^i^^i'^^i^°i^'^i^^ as a 
mesoscopic simulation technique (not only as a thermostat) to simulate hydrodynamic be- 
havior as well as the rheological properties of complex fluids when thermal fluctuations are 
important. The fluid is modeled in terms of mesoscopic particles known as dissipative par- 
ticles that are large quasi-particles, which evolve in the same way as MD particles do, but 
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with different inter-particle interactions that allow for much longer time steps. The forces 
between each pair of dissipative particles is made up of a conservative force, a dissipative 
force and a random force similar to the one given by Eq. but each of which is pairwise 
additive. These forces in general conserve total momentum and have a spatial range given 
by the cut-off distance r^. In this case, the dissipative and random forces are the main 
ingredients for hydrodynamic interaction and are commonly given by, 

Ffj = -m^uj^{rij){rij ■ Vij)rij, F^j = mau^{rij)Cijrij, (13) 

where 7 and a are constants, uj^{r) and uj^{r) are weight functions that vanishes for r > r^, 
Vij and Vij are relative velocity and distance between the pairs, respectively, Qj is a random 
noise term with zero mean. The random number can be sampled from either a Gaussian 
or a uniform distribution, in each case with a variance of unity. Espafiol and Waven^ 
have showed that although the weight function can be chosen arbitrarily, it must satisfy the 
relation uj^{r) = [uj^{r)]'^. The friction coefficient 7 and the noise amplitude a are related 
by the fluctuation/dissipation theorem, ma^ = 2'~fkBT. 

In the present study, however, we used DPD to thermostat the system to take advantage of 
the fact that it conserves the hydrodynamic interaction between the particles. The function 
used for uj^{r) and uj^{r) in this study have the following forn>^ 

uj^{r) = [cu (r)]^ = < (14) 
I 0, r > 

where we used the same cutoff as for the Lennard- Jones interaction cutoff between 
monomers, r'^ = = 2.5a. 

All the simulations were run using the massively parallel code LAMMPS.— The equations 
of motion were integrated with a velocity- Verlet algorithm with a time step At = 0.012r for 
the interdiffusion study and At = 0.009r for the bulk equilibrium measurement of the self 
and corrected diffusion constants, where r = m(o"/e)^/^. The smaller time step was used to 
assure that the system was stable even for large 7 (see Fig. |2I). All the simulations were car- 
ried out at a temperature of T = e/ks and pressure P ^ 0. To determine the fugacity f, the 
LADERA grand canonical molecular dynamics (GCMD) code was used.— During the course 
of an equilibrium molecular dynamics simulation at the appropriate solvent concentration, 
the energy of inserting a solvent particle at random locations was sampled. 
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The simulated system for the interdiffusion study consists of a rectangular cell, which 
is periodic in x and y direction but not in z as shown in Fig. ^ This initial configuration 
was generated in two steps. First, a polymer melt system {L^ = Ly = 60cr and thickness 
Lz — 9Qa) was equilibrated between two walls at pressure P ~ 0. Then the top of the 
box was extended and the solvent molecules were placed in contact with the polymer melt. 
Furthermore to keep the pressure of the system constant, a small vapor phase is added 
between the top wall and the solvent particles. The polymer melts in this study consisted 
of chains of length N = 500 or 50 monomers. The total number of polymers monomers in 
all cases was 300,000. The solvent consisted of either 230,000 monomers or 193,000 dimers. 

For studying self and corrected diffusion as a function of solvent concentration, the system 
consists of an equilibrated polymer solvent mixture in a cubic cell, which is periodic in all 
three directions. The total number of polymer monomers Np = 50, 000 in 100 chains each 
with 500 monomers. The number of solvent monomers Ng in the system was varied from 500 
monomers (dilute case) to 150,000 monomers. A pure solvent system of 50,000 monomers 
was also simulated. As the concentration of solvent is varied, the pressure in the system 
is kept constant by allowing the system relax to the corresponding volume, thus fixing the 
volume for the measurement of the diffusion constant. Only a few hundred thousand MD 
timesteps are required to calculate the self diffusion constant Ds{c) though more than ten 
times this is required for the corrected diffusion constant Dc{c). The fugacity calculation 
requires approximately a million MC cycles and a few hundred thousand MD timesteps 
especially at low solvent concentration. 

III. DIFFUSION COEFFICIENTS 

A. Dependence of self diffusion on the strength of the dissipative force 7 

A set of simulations were performed to determine the dependence of the diffusion on 7. 
The self diffusion constant, Ds{c), of the solvent in the system is calculated using Eq. 01 
The results for -Ds(l) are shown in Fig. El for the two thermostats. 

We find that the values of -Ds(l) from the two thermostats strongly depend on 7 and 
both thermostats show similar dependence. However, for a given value of 7, -Ds(l) from 
the two thermostats are not necessarily expected to be equal and the agreement for the two 
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FIG. 2: Dependence of self diffusion constant -Ds(l) for the pure solvent on the damping coefficient 
7 for the Langevin (a) and DPD with = 2.5(T (o) and = 2. Oct (■) thermostats for T = e/A;^ 
and P ~ 0. The solid lines are a guide to the eye. Error bars are ib0.002(T^/T. 

thermostats for = 2.5o" is coincidental. As shown in Fig. El -0^(1) depends not only on 7 
but also on the cutoff in Eq.^J For small values of 7, -Ds(l) decreases exponentially, while 
for large values of 7, -Ds(l) decreases as 7"^ as expected since the random and dissipative 
forces dominate the interaction between particles and the diffusion becomes Brownian. Note 
that for large values of 7 the Langevin thermostat becomes unstable for the value of timestep 
used, At = 0.009r. The strong dependence of -Ds(l) on 7 is often ignored specially for DPD 
particles where large values of 7 and large time steps are usually taken. In the rest of our 
simulation we choose 7 = O.lr^^ so that the interdiffusion is only weakly affected by the 
random and dissipative forces. 

B. Concentration dependence of diffusion constants 

The self diffusion, -Ds(c), and corrected diffusion, -Dc(c), constants are calculated from 
Eq. 01 and d respectively. The calculated values as a function of solvent concentration 
are shown in Fig. 01 For comparison, the two thermostats with 7 = O.lr^^ were used to 
compute the self diffusion. The corrected diffusion constant was determined using the DPD 
thermostat. The self diffusion values from the two thermostats are the same within the 
error of the simulation. In general, for low solvent concentration both the self and corrected 
diffusion constants show weak dependence on concentration. The self diffusion constant then 
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FIG. 3: Dependence of diffusion constants on solvent concentration for Ds{c) using the Langevin 
thermostat (a) and DPD thermostat (o), and Dc{c) (□) using the DPD thermostat and D{c) (■) 
from Darken equation Eq. ^ For c < 0.09, points for Dc{c) and D{c) overlap. The solid lines are 
a guide to the eye and all results are for polymer chain length N = 500 and ei2 = e. Error bars 
are ±Qm2a'^/T. 

increases slowly for intermediate solvent concentration, but a sharp increase is observed for 
larger concentration. 
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FIG. 4: Fugacity as a function of solvent concentration for monomers in a polymer melt of chain 
length N=500 and ei2 = e. Error bars are ±0.004. 

To calculate the diffusivity D{c) from the Darken equation, Eq. 0J the thermodynamic 
factor (91n//(91nc is also required. The fugacity of the solvent was calculated using the 
GCMD simulation method^ and is shown in Fig. |3] Numerical differentiation of this data 
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gives the required thermodynamic factor as a function of solvent concentration. The ther- 
modynamic factor decreases monotonically with concentration opposite to that obtained for 
solvent in a zeolite or other porous material due to the fact that the solvent swells the poly- 
mer, making it easier not harder to insert solvent monomer as the solvent density increases. 
The diffusivity, D{c), calculated from the Darken equation is shown in Fig. |21 Note that 
D{c) and -Ds(c) are equal at the dilute limit and show a similar dependence on concentration, 
though -D(c) increases slower with increasing concentration than -Ds(c). 

IV. INTERDIFFUSION 

The initial setup for the interdiffusion studies of a solvent into an equilibrated polymer 
melt is shown in Fig.^ The density profile, as usually used in experiments, of both polymer, 
Pp and solvent, Ps as a function of time for the two thermostats is shown in Fig. Elfor a 
monomer solvent diffusing into a polymer melt of chain length = 500. As the solvent 
diffuses into the polymer, the polymer relaxes and the boundary is smeared out. The density 
profile at different times for the two thermostats have the same shape differing slightly only 
in time scale. The rate at which the solvent penetrates into the polymer can be determined 
either taking a particular value of Ps to define the depth of penetration or by the weight gain 
by the polymer system as a function of time. The weight gain by the polymer versus t^^"^ is 
shown in Figure El Both thermostats give the same result confirming that the penetration 
increase with t^^"^ in agreement with the Fickian diffusion. This can be further justified by 
plotting the solvent density profiles of Fig. El for the DPD thermostat function of 

as shown in Fig. [71 The profiles collapse on to a single master curve confirming the 
Fickian behavior. 

As a first approximation we treat D{ clS cL constant Dq as often done experimentally 
and fit the solvent concentration profiles of Fig. (3 to the erf function given by Eq. |21 As 
seen from in Fig. |7|^a) (solid line), the erf function fit the concentration profile reasonably 
well particularly in the low solvent concentration, though overestimates slightly the rate the 
monomer solvent enters the polymer film near the interface. The diffusivity extracted from 
the fit in Fig. IHa) is Dq ~ 0.033 ± 0.002 a^/r independent of time. This value is within 
error bars in agreement with D{c) in the dilute limit, -D(O) = 0.030 ± 0.002 a^/r from Sec. 
III. 
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FIG. 5: Solvent ps and polymer pp concentration profiles as a function of time starting from t = 
and plotted every 2400 r. The solvent is diffusing into the polymer from the right side and (a) is 
for the Langevin thermostat and (b) is for the DPD thermostat. 




FIG. 6: Weight gain for three solvent-polymer systems, monomer solvent for polymer chain length 
N=500 (a) and N=50 (□) and dimer solvent for chain length N=500 (o) using DPD thermostat. 
Similar results for N=500 and monomer solvent were obtained for the Langevin thermostat. 

Using the predicted form of D{c) shown in Fig. El from the Darken equation, Eq. ^can 
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FIG. 7: Solvent concentration profiles are plotted as a function of the scaling variable zt^^/^, (a) 
for a monomer solvent that was shown in Fig. [S] and (b) for a dimer solvent. The solid line in both 
cases is a theoretical fit using Eq. |21for the profiles. The fit using D{c) from Fig. El is shown as 
dashed line. 

be solved numerically. D{c) can be fit to 0.031/(1 — 1.15c). We find the result shown by the 
dashed line in Fig.[7|^a), which fits the simulation result very well at low solvent concentration 
and overestimates slightly more than constant diffusivity at large concentrations. This 
deviation at large concentration can be partly attributed to the swelling of the polymer 
which is neglected in Eq. ^ In Fig- IZ(b) we show solvent concentration profiles for the 
case of dimer solvent and due to the smaller density difference between the solvent and the 
polymer melt there is less swelling. Since it is difficult to calculate the fugacity for the dimer 
solvent case we did not calculate the diffusivity as a function of concentration. In this case 
a constant diffusivity fit (solid line) gives very good agreement with the simulation results 
even at large concentration. The diffusivity extracted from this fit is Dq = 0.017 ± 0.001 
(T^/r, which also equaled -Ds(O). Note that the increase in Ds{c) from the dilute to the pure 
solvent limit is about a factor of 6 for the monomer solvent case and a factor of 3.5 for the 
dimer solvent case. An extrapolation of the fit to the monomer D{c) gives a factor of 2.9 
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from the dilute to the pure solvent limit which suggest that for the dimer case D{c) should 
increase very slowly with concentration, which further supports the assumption of constant 
diffuivity. 

Our result agrees with the recent hypothesis^^^ based on two-dimensional lattice gas 
automation simulation, that the nature of the diffusion is not related to the solvent con- 
centration gradient within the system but to the diffusivity gradient {dD{c)/dc), and the 
standard Fickian diffusion only occurs when dD{c)/dc ~ 0. This justifies the experimental 
fits to the erf function. 

We also explored interdiffusion in a number of different ways, including varying ei2 
(polymer-solvent interaction parameter) and polymer chain length. First consider the effect 
of varying €12- For 612 = 0.8e the solvent did not diffuse into the polymer melt for the time 
scale of our simulation. In this case the solvent can be considered to be a poor solvent. On 
the other hand, for €12 = 1.2e little change from the previous Fickian diffusion behavior was 
observed. For this case the diffusion constant extracted from the error function fit is slightly 
larger, Dq = 0.037 ± 0.003 a^/T,than for = 1.0. 

To study the effect of the chain length of the polymer we decreased the chain length from 

= 500 to 50. As seen from Fig. IHlthe diffusion process did not change as expected. The 
diffusion constant Dq for a monomer solvent extracted from the error function fit for N =50 
is Do = 0.035 ± 0.002 a^/r compared to D,(0) = 0.032 ± 0.002 a^/r. 

V. SUMMARY 

In this work large scale molecular dynamics and grand canonical Monte Carlo simulation 
techniques were used to study the interdiffusion of a solvent into a polymer melt. The self 
and corrected diffusion constants as a function of solvent concentration were determined 
separately and compared to those obtained from the interdiffusion studies using the Darken 
equation. For low and intermediate solvent concentration, both Ds{c) and -Dc(c) increased 
slowly with solvent concentration and were equal within the error of the simulation. For 
larger concentrations both diffusion constants increased rapidly with the corrected diffusion 
constant increasing significantly faster than the self diffusion constant. Because the solvent 
swells the polymer, the thermodynamic factor din f /dine decreased with increased solvent 
concentration resulting in a diffusivity D{c) which is essentially constant for low and inter- 
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mediate concentration and increased less rapidly at high concentration than both Dc{c) and 
Ds{c). The observed dependence of D{c) with concentration is opposite to that of a zeolite 
where D(c) increases more rapidly with concentration than Dc{c). 

Fickian diffusion behavior was observed for solvent absorption into polymer melt for all 
cases studied. This was verified by the t^^"^ dependence of the weight gain by the polymer 
system and thus the diffusion process can be considered to be Fickian. The concentration 
profile of the solvent fit an error function derived from Fick's second law for constant diffu- 
sivity. The diffusivity found from this fit was found to be independent of time and is equal 
to the self diffusion constant Ds{0) at the dilute limit. Even though D(c) is not constant 
over the entire range of concentration, since it varied little in the low concentration region 
relevant for interdiffusion, assuming D{c) constant is a very good approximation. 

We also studied the dependence of the interdiffusion on the polymer-solvent interaction 
strength, the chain length of the polymer and the chain length of the solvent. When the 
interaction parameter was slightly lowered from the neutral case the solvent did not diffuse 
into the polymer on the time scale of our simulation. On the other hand, increasing the 
polymer-solvent interaction parameter by the same amount did not considerably affect the 
diffusion process. The diffusion process was not also affected by a change in the chain length 
of the polymer. 

Future work will study the crossover from Fickian to non-Fickian diffusion as the state 
of the polymer changes from a melt to a glass. 

VI. ACKNOWLEDGMENTS 

We are thankful to A. P. Thompson for providing the GCMD code and for helpful dis- 
cussions. Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed 
Martin company, for the United States Department of Energy's National Nuclear Security 
Administration under Contract No. DE-AC04-94AL85000. 



1 J. S. Vrentas, J. L. Duda, and W. J. Huang, Macromolecules 19, 1718 (1986). 
^ J. Sonnenburg, J. Gao, and J. H. Weiner, Macromolecules 23, 4653 (1990). 
3 H. Takeuchi, R.-J. Roe, and J. E. Mark, J. Chem. Phys. 93, 9042 (1990). 

15 



^ p. V. K. Pant and R. H. Boyd, Macromolecules 25, 494 (1992). 

^ F. Miiller-Plathe, S. C. Rogers, and W. F. van Gunsteren, Macromolecules 25, 6722 (1992). 

^ R. M. Sok, H. J. C. Berendsen, and W. F. van Gunsteren, J. Chem. Phys. 96, 4699 (1992). 

A. A. Guscv, S. Arizzi, U. W. Sutcr, and D. J. Moll, J. Chem. Phys. 99, 2221 (1993). 

^ Y. Tamai, H. Tanaka, and K. Nakanishi, Macromolecules 27, 4498 (1994). 

9 R. H. Gee and R. H. Boyd, Polymer 36, 1435 (1995). 

1° D. Hofman, J. Ulbrich, L. Fritz, and D. Paul, Polymer 37, 4773 (1996). 

11 T. Li, D. O. Kildsig, and K. Park, J. Controlled Release 48, 57 (1997). 

12 F. Miiller-Plathc, J. Chem. Phys. 108, 8252 (1998). 

13 M. C. Griffiths, J. Strauch, M. J. Monteiro, and R. G. Gilbert, Macromolecules 31, 7835 (1998). 

14 O. Hahn, D. A. Mooney, F. Miiller-Plathe, and K. Kremer, J. Chem. Phys. Ill, 6061 (1999). 
1^ N. F. A. van der Vegt, Macromolecules 33, 3153 (2000). 

16 M. L. Greenfield and D. N. Theodorou, Macrol. 34, 8541 (2001). 

1^ E. Tocci, D. Hofmann, D. Paul, N. Russo, and E. Drioh, Polymer 42, 521 (2001). 

18 S. Y. Lim, T. T. Tsotsis, and M. Sahimi, J. Chem. Phys. 119, 496 (2003). 

1^ C. J. Burning, M. M. Hassan, H. M. Tong, and K. W. Lee, Macromolecules 28, 4234 (1995). 

20 M. M. Hassan and C. J. Burning, J. Polym. Sci. 37, 3159 (1999). 

21 M. Sanopoulous and F. StamatiaHs, Polymer 42, 1429 (2001). 

22 D. F. Stamatialis, M. Sanopoulous, and J. H. Petropoulos, Macromolecules 35, 1021 (2002). 

23 M. Tsige and P. L. Taylor, in preparation (2003). 

24 W. Jilge, I. Carmesin, K. Kremer, and K. Binder, Macromolecules 23, 5001 (1990). 

25 H. P. Beutsch and K. Binder, J. Chem. Phys. 94, 2294 (1991). 

26 J. Crank, The Mathematics of Diffusion (Oxford University Press, Oxford, 1975). 
2"^ N. J. M. Kuipers and A. A. C. M. Beenackers, Chem. Eng. Science 48, 2957 (1993). 
28 N. L. Thomas and A. H. Windle, Polymer 22, 627 (1981). 

2^ I. Quijada-Garrido, M. F. de Velasco-Ruiz, and J. M. Barrales-Rienda, Macromol. Chem. Phys. 
201, 375 (2000). 

30 N. L. Thomas and A. H. Windle, Polymer 21, 613 (1980). 

31 N. L. Thomas and A. H. Windle, Polymer 23, 529 (1982). 

32 C.-Y. Hui and K.-C. Wu, J. Appl. Phys. 61, 5129 (1987). 

33 L. S. Barken, Trans. AIME 175, 184 (1948). 

16 



R. Sharma and K. Tankeshwar, J. Chem. Phys. 108, 2601 (1998). 
35 E. J. Maginn, A. T. Bell, and D. N. Theodorou, J. Phys. Chem. 97, 4173 (1993). 
3^ D. N. Theodorou, Diffusion in Polymers (Marcel Dekker Inc., New York, 1996). 
3^ A. I. Skoulidas and D. S. Sholl, J. Phys. Chem. B 105, 3151 (2001). 

38 M. Chandross, E. B. W. Ill, G. S. Crest, M. G. Martin, A. P. Thompson, and M. W. Roth, J. 
Phy. Chem. B 105, 5700 (2001). 

39 J. M. D. MacElroy and K. Raghavan, J. Chem. Phys. 93, 2068 (1990). 

^° A. P. Thompson and G. S. HefFelfinger, J. Chem. Phys. 110, 10693 (1999). 

K. Krcmcr and G. S. Grest, J. Chem. Phys. 92, 5057 (1990). 
^2 G. S. Grest and K. Kremer, Phys. Rev. A 33, 3628 (1986). 
^3 B. Diinweg, to be published (2003). 

T. Soddemann, B. Diinweg, and K. Kremer, Phys. Rev. E to be published (2003). 
^5 P. J. Hoogerbrugge and J. M. V. A. Koleman, Europhys. Lett. 19, 155 (1992). 
^6 P. Espanol and P. Warren, Europhys. Lett. 30, 191 (1995). 
^7 R. D. Groot and P. B. Warren, J. Chem. Phys. 107, 4423 (1997). 

48 J. B. Gibson, K. Chen, and S. Chynoweth, J. Colloid and Int. Sci. 206, 464 (1998). 

49 J. B. Gibson, K. Zhang, K. Chen, S. Chynoweth, and C. W. Manke, Molecular Simulations 23, 
1 (1999). 

50 P. Malfreyt and D. J. Tildesley, Langmuir 16, 4732 (2000). 

51 J. A. Elliott and A. H. Windle, J. Chem. Phys. 113, 10367 (2000). 

52 M. Ripoll, M. H. Ernst, and P. Espanol, J. Chem. Phys. 115, 7271 (2001). 

53 S. J. Plimpton, J. Comput. Phys. 117, 1 (1995). 

54 M. Kuntz and P. Lavallee, J. Phys. D: Appl. Phys. 36, 1135 (2003). 



17 



